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ABSTRACT 

The detection of pulsations from an X-ray binary is an unambiguous signature of the 
presence of a neutron star in the system. When the pulsations are missed in the radio band, 
their detection at other wavelengths, like X-ray or gamma-rays, requires orbital demodulation, 
since the length of the observations are often comparable to, or longer than the system orbital 
period. The detailed knowledge of the orbital parameters of binary systems plays a crucial 
role in the detection of the spin period of pulsars, since any uncertainty in their determination 
translates into a loss in the coherence of the signal during the demodulation process. In this 
paper, we present an analytical study aimed at unveiling how the uncertainties in the orbital 
parameters might impact on periodicity searches. We find a correlation between the power 
of the signal in the demodulated arrival time series and the uncertainty in each of the orbital 
parameters. This correlation is also a function of the pulsar frequency. We test our analytical 
results with numerical simulations, finding good agreement between them. Finally, we apply 
our study to the cases of LS 5039 and LS 1 +61 303 and consider the current level of uncer- 
tainties in the orbital parameters of these systems and their impact on a possible detection of 
a hosted pulsar. We also discuss the possible appearance of a sideband ambiguity in real data. 
The latter can occur when, due to the use of uncertain orbital parameters, the power of a puta- 
tive pulsar is distributed in frequencies lying nearby the pulsar period. Even if the appearance 
of a sideband is already a signature of a pulsar component, it may introduce an ambiguity in 
the determination of its period. We present here a method to solve the sideband issue. 

Key words: stars: neutron, pulsars, gamma-rays: observations 



1 INTRODUCTION 

A few High Mass X-ray Binaries (HMXBs) ha ve been detected to 
emit GeV and TeV photons. LS I +61 30 3 (seelAbdo et alj|2009bl . 
lAlbert et ail 120091 TAcciari et af 
5039 (seelAbdo et all 



2009c 



63/LS 2883 (see lAbdo et al 
J0632+057 



(see 



2oTil 



Albert et al. 2006), LS 



Aharonian et al 



2011a 



Aharonian et al 



20061) PSR B 1259- 



Aharonian etal]|2005l) , HES S 



20071 iBongiorno etah 2011) 



and th e most recently discovered 1FGL J 101 8. 6-5856 lAbdo et al.l 
d2012l) are all examples of these systems. These former systems are 
generally referred to as gamma-ray binaries. Other objects like Cyg 
X-l and Cyg X-3 have been observed to flare in gamma-rays, but 
thei r emission is neither dominant nor pe rsistent at these energies 
(see lAlbert et al] 20071 lAbdo et alj2009ah . 

These HMXB systems are composed by a massive OB 
or Be star and a compact object, the nature of which is 
in general unknown. In the case of Cyg X-l (and less se- 
curely of Cyg X-3) it is likely, however, that their compact 
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objects are black holes surrounded by an accretion disk filled 
by matter captured from the massive star. On the other hand, 
among the gamma-ray emitting HMXBs, the compact object 
of PSR B1259 -63/LS 2883 i s a rad io puls ar with a spin pe- 
riod o f 48 ms I Johnston et all i 19991). After [ Maraschi & Treves! 



19811). e.g.. 



Dubusl d2006h. Isieroowska & Torres! J2008l) . and 
IZdziarski. Neronov. & Chernvakov al d2010h . among others, pro- 
posed detailed theoretical models to explain the emission of these 
gamma-ray binaries as due to the interaction of the relativistic par- 
ticle wind from the pulsar with the wind of the massive star, or via 
processes in the pulsar wind zone directly. Recently, following the 
detection of a very shor t, magnetar-like bur st coming from the di- 
rection of LS I +61 303. lTorresetaI] l l2012l) developed on a model 
based on assuming the existence of a high magnetic field - long pe- 
riod pulsar in this system. The dichotomy on the nature of gamma- 
ray binaries is a tre nding topic of high-energy astrophysics (see, 
e.g., [Mirabel 2006). Obviously, the detection of pulsations from 
them would lead to a clear and unequivocal solution. 

Deep searches for pulsations in radio frequencies have been 
performed especially for LS I +61 303 and LS 5039 (e.g. 
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iMcSwain etal]201lh . without success. The lack of radio pulsation 
can be explained by the dense environment of the massive star. In- 
deed, for PSR B 1259-63/LS 2883 (the largest system, with a period 
of ~ 4 years) the radio pulsations disappear at t he periastron dis- 
tance l ljohnston et alj[l99^ . I Johnston et alj|2005l) . Since the orbits 
of LS I +61 303 and LS 5039 are much smaller than that of PSR 
B 1259-63, it is natural to expect the radio emission of their compact 
object being always affected by free-free absorption and disper- 
sion. Searches for pulsations were performed also with X-ray data 
for LS 1+61 303, and LS 5039 and the deepest upper limits on the 
pulse fr actio n were put in the w orks bv lRea et al.l j2010h .l ReaetaU 
l l2011al) . and lRea et Note again that for the only firmly 

established TeV binary containing a p ulsar, PSR B1259-63, X - 
ray pulsations were not detected either JChernvakova et alj|2009h . 
pointing to an X-ray emission being dominated by wind-wind or 
intra-wind shocks. 

The Fermi satellite, and its main instrument on board, the 
Large Area Telescope, is continuously surv eying the sky since its 
launch in June 2008 dAtwood et alJ.l2009l) . This experiment of- 
fers a good opportunity to perform searches for gamma-ray pul- 
sations from binaries. However, given the dim character of gamma- 
ray fluxes, and the paucity of counts, large integration times are 
needed for pulsation searches. The event arrival times need orbital 
demodulation. An uncertain knowledge of the orbital parameters, 
in this gamma-ray case, or at any other frequencies or systems 
where demodulation is needed, would lead to lose the coherence 
of the pulsed signal. The most accurate measurements of the or- 
bital parameters for LS 5039 and LS 1+61 303 are derived fitting 
the Dop pler shift of optical spectral lines emitted by the massive 
star ( s ee|Casares et al.l2005l ICasares et alj|2005l iGrundstrom et al.l 
12001 lAragona et alj|2009ir This technique led so far to uncertain- 
ties of the order of 1% to 10%. Inverting the problem, one may 
ask how well should a certain orbital parameter be known in order 
to secure that a pulsed signal is not lost through the demodulation 
process. 

The purpose of this work is to analytically study how much the 
uncertainties on the orbital parameters affect the power spectrum of 
a putative pulsed signal from an X-ray binary. 

The paper is organized as follows. Section 2 presents our ap- 
proach to the problem, using perturbation theory. Section 3 intro- 
duces the perturbation functions, and a new set of variables that 
simplify the treatment of the problem. In Section 4 and 5, we 
compute the probability density function of the perturbed emission 
times, and for each orbital parameter, the impact of their uncer- 
tainty in the power spectrum. Section 6 provides constraints on the 
uncertainty over each of the orbital parameters such that the loss 
in the power of the pulsar signal is smaller than a given value. Fi- 
nally, the discussion applies our results to LS 5039 and LS I 61 303, 
together with providing numerical simulations that validate our re- 
sults and conclusions. Numerical validation of our conclusions is 
given via simulations and subsequent timing analysis of pulsed sig- 
nals from different binaries, where the knowledge of the orbital 
parameters is blurred ad-hoc. 



2 THE PROBLEM OF PULSE EXTRACTION 

2.1 Numerical problems for a blind search approach 

We start showing how a blind search approach for pulsars in binary 
systems is generally doomed. 

For a blind search of isolated pulsars, a FFT should be per- 
formed with a number of frequency bins (Nf) equal to Xobs-pMAx, 



where T b s is the viewing period, and Fmax the highest frequency 
searched. Many trials are needed to correct for the first derivative 
of the frequency, Fl. The step size in Fl should not be larger than 
1 /Tote m order to keep the signal power within a single frequency 
bin. The number of Fl trials, Nfi (that is proportional to T^ bs ) 
would be enormous for viewing periods lasting few years. Then, 
even for isolated pulsars, performing a blind search is computa- 
tional demanding. To face this issue different techniques were pro- 
posed, the most suc cessful so far bein g the time-differencin g of 
lAtwood, et ail ( 120061) . and the method of lPletsch. et"al1 j2012h . 

In the case of a pulsar in a binary system, many trials would be 
needed also for each of the orbital parameters. We could arbitrary 
choose to cover the uncertainty ranges of the orbital parameters (we 
consider just 5 parameters) with an equal number of trials N p each. 
The total number of trials would be Nt = Np x Nfi- 

To have an idea on the n umbers involved, we can c onsider a 
blind search similar to that in lSaz Parkinson et al.l ( I20101) . but ap- 
plied to binary systems. In that case, pulsations are searched in one 
year of Fermi-LKT data using the time-differencing technique with 
time windows of 6 days, up to 64 Hz frequency. The number of 
frequency bins in the FFT were Nf = 2 26 , and the trials in the 
frequency derivative were Nfi = 2000. In the case of binary sys- 
tems, an hypothesis of just 10 trials for each orbital parameter leads 
to a total number of trials equal to Nt — 2 x 10 s . 

And yet, is N p = 10 a guarantee that the pulsar detection does 
not get lost, regardless of its frequency? How can we understand 
whether the uncertainty range is oversampled or not by the choice 
of the number of trials? Is it possible to optimize the trials in order 
to minimize the required CPU time to run the analysis? It is not 
possible to answer all these questions without studying how the 
parameter uncertainties affect the results of the periodicity search. 

2.2 The analytical approach 

The photons emitted by a pulsar in a binary system experience sev- 
eral delays travelling towards the observer. These delays are: the 
dispersive delay due to the Interstellar Medium (A/s), the prop- 
agation and relativistic delays within the Solar System (A©), and 
the corresponding delays accounting for the geometry in the bi- 
nary system itself (As). The so-called timing formula correlates 
the photon emission time in the pulsar reference frame (if sr ) with 
the photon arrival time to the observer (t° bs ) as 

tT =t°« s -Ao-Ajs-Ab. (1) 

The delay introduced by the motion of the pulsar around its com- 
panion star (As) is mainly due to the Romer delay (Ar), 

A B ~A R = 

A [sinW (cos_B - e) + y/l - e 2 ■ cosVK sinE 1 ] , (2) 

where A is the projection of the semi-major axis on the plane per- 
pendicular to the observer's line-of-sight, W is the longitude of 
periastron, e is the eccentricity, and E is the eccentric anom aly 
(see for example iBlandford and Teukolskvl ll976.Eqn2.2l . or 
ICamenzindll2007. chapter 6.5.2. pag. 259l) . Post-Newtonian effects 
will be neglected in this work. A schematic view of an orbit and 
its parameters is shown in Figure Q] while Table Q] lists the main 
variables. 

Assuming that a pulsar is part of a binary system, the detec- 
tion of periodic signals from it can be achieved by calculating the 
power spectrum of the emission time series in Eq. (1). In practice, 
however, one is faced with the problem that the uncertainties in the 
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line of nodes 



Y 

Figure 1. Schematic view of an orbit and its parameters. The gray plane is 
the reference plane perpendicular to the observer line of sight, and i is the 
inclination of the plane of the orbit respect to it. a and b are the semi-major 
and the semi-minor axis, respectively. F is the focus of the orbit that hosts 
the baricenter of the binary system, while C is its center. The distance CF is 
equal to the semi-major axis times the eccentricity (a ■ e). W is the longitude 
of the periastron, and v is the true anomaly. 



estimation of the delays lead to a wrong calculation of the emis- 
sion times (tew), affecting the results of the analysis. In this work, 
we focus on studying the impact of the uncertainties in the orbital 
parameters on the ability for detecting pulsations. For this reason, 
it will be assumed that the delays due to the interstellar medium 
(Ajs), and to the propagation in the solar system (A@) are exactly 
known (or rather, that their uncertainty is much smaller in compari- 
son to Ah). The Romer delay in the binary system (Eq. 2), instead, 
is considered wrongly estimated (Ar w ). 

Comparing the timing formula (Eq. 1) in the cases of a cor- 
rect and a wrong calculation of the Romer delay (An and Ar w 
respectively), the relation between the correct (t e ) and wrong (tew) 
emission times turns out to be 



t, 



SA R , 



(3) 



where SAr = Ar w — Ar. A mathematical interpretation of 
Eq. (O is simply that the times t ew are the correct emission times 
perturbed by SAr. At first order approximation, the impact of the 
errors in the orbital parameters on the Romer delay calculation, can 
then be taken into account by writing the perturbation factor (SAr) 
as 



(4) 



where the sum is over the orbital parameters p = 
{A, W, e, P or b, To}, and To is the epoch of the periastron. 
Eq. © allows us to factorize the problem, and to evaluate how the 
timing analysis is affected by each orbital parameter separately 
from the others. With this aim we introduce the perturbation 
functions, f p , relative to each of the orbital parameter, p, as 



U = -^dp. 



(5) 



3 THE PERTURBATION FUNCTIONS 

In Eq. lO the perturbation functions are defined as partial deriva- 
tives of the Romer delay (A_r), which is in turn expressed in 



Eq. {2]l. It is possible to rewrite this latter formula in a more com- 
pact form, using the following trigonometric identity: 



a ■ cos x + b ■ sin x = c ■ sin(a; + (f>) , 
<f> — atan(a/fc), 



(6) 



C= \/a 2 + b' 2 . 

Applying it to cos E and sin E in Eq. J2}, the Romer delay is then: 

A R = M(A, e, W) ■ sin(P + 0(e, W)) - Q(A, e, W), (7) 
where 



M = A*J\ - e 2 cosW, 



<j> = atan(tan(W / ")/\/l-e 2 )), 
Q = e A sinW. 



(8) 



To evaluate the perturbation functions, we calculate the partial 
derivatives of Eq. Q for a generic orbital parameter p, 



dAt 



dM 



dp dp 



sin(T - 



,'dE d4>\ dQ 
+ M.cos(T + ^).(- + ^j-^. 



(9) 



Here again, we can apply the trigonometric identity in Eq. ((6) to 

sin(£' + <f>) and cos(£' + (f>), obtaining: 

8A R 



dp 



K ■ sin(S + <t> + ip) - C, 



(10) 



from which we can define 

5k 



Kdp 



{d E + 8J, 
dp dp 



+ 



dM 



dp 



dp. 



at an 



(ii) 



Sc — Cdp = ^Q-dp. 

dp 

Finally, a compact formula for the perturbation function is 



/ = 



8 A f 
dp 



■dp = 8k ■ sin(i5 + (j> + tp) — Sc. 



(12) 



Since this formula is valid for all the orbital parameters, we omitted 
the sub-index p in Eq. J 1 2b . From a mathematical point of view, we 
have done a transformation from the canonical orbital parameters 
{A, W, e, Porb, ?o} and their errors {dA, dW, de, dP olb , T }, to 
the parameters {8k, <f>, ip, 5c}. 

The behavior of the perturbation function, at each instant, is 
given by its dependence on the eccentric anomaly E, and it is con- 
nected to the emission times t e by the relation 



E-e-smE = n oTh (t e -T ), 



(13) 



where Q rb = 27r/P or b is the system frequency. This equation 
provides the emission times as a function of E. The inverse func- 
tion (E — E(t e )) can not be expressed analytically using Eq. J 13b . 
However, it can be approximated, and in order to express the pertur- 
bation function dependence with t e and do further analytical steps, 
we will simply assume that E ~ Qorbt e - This implies, 



f(t e ) ~ 5k ■ sin(n orb te + (f> + i>) — 8c. 



(14) 



Appendix [B] gives an assessment of the approximation made to 
reach the latter formula. Note that 8k and ip can also be functions of 
t e , because their defining formulae (see Eq.ll 111 contain the partial 
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Table 1. Meaning of the variables used in the paper 



Variable Meaning 



,psr . 

t% or t e 


Photon emission time in the pulsar reference frame 


xobs 
L a 




Photon arrival time to the observer 


-* 

veu) 




Perturbed emission time: photon emission time in the pulsar reference frame affected by the errors on the orbital parameters 


As 




Romer delay 






Evaluation of the Romer delay affected by the errors on the orbital parameters 


A 




Projection of the semi-major axis on the plane perpendicular to the observer line of sight 


W 




Longitude of periastron 


e 




Eccentricity 


Porb 




Orbital period 


^orb 


= 27r/P orb 


Angular frequency of the binary system 


E 




Eccentric anomaly 


To 




Epoch of the periastron 


f p or 


/ 


Perturbation function, defined in Eq. j5J, and described in Section 3 


{6k, < 


i>, 4>, 8c} 


Set of transformed orbital parameters. <f> is defined in Eq. (8), the other three in Eq. Jilt 


rp 

J obs 




Duration of the observation 


OJQ = 




Angular pulsar frequency 



derivatives dE/dp. By using Eq. d!3t . we can see that the partial 
derivative is null when the orbital parameter p incarnates into A, 
or W. In these two cases, 8k and ip are constants. In contrast, the 
partial derivative dE/dp is a function of E when p is either e or 
Port,, since 8k and tj) are function of E, and consequently, also of 
t e . This difference in the dependence of 8k and ip will lead to a 
different treatment of the problem, as described below in Section 6. 



4 THE PROBABILITY DENSITY FUNCTION OF THE 
PERTURBED EMISSION TIMES 

The expectation value of the power spectrum can be evaluated us- 
ing the probability density function (pdf) of the phases assigned to 
each photon (see appendix lAl). In order to evaluate the power spec- 
trum of the perturbed emission times t ew , we shall calculate their 
pdf (Pt BW ), as well as that of the phases assigned to them {Pg). 

Since the times t ew are correlated with the correct emission 
times (t e ) by Eq. l[3}, which we rewrite here as 

t ew =t e - f(t e ), (15) 

Pt ew can be calculated if Pt c is known. Appendixlclshows that for 
all realistically observable binary systems, t ew is a monotonic in- 
creasing function of t e . This makes the calculation of its pdf easier. 
Indeed, its pdf is: 

P.M = uJ^L y (i6) 

where here (and hereafter) U indicates a normalization factor, and 
a / represents the first derivative respect to the emission time t e . 
Finally, to compute the power spectrum of the perturbed time series 
t em at a frequency uj one has to calculate the phases, defined as 

9 = wW (17) 

Since to acts like a constant in Eq. J17t , the pdf of the phases has 
the same form of Eq. J16b . i.e., 

Mt * ] = U ^T§7) (18) 

The fundamental features of the power spectrum of the per- 
turbed time series t ew can be derived assuming that the signal emit- 
ted by the pulsar in its reference frame is sinusoidal with frequency 



wo, 

Pt K (te) = 1 + sin(oj ie). (19) 
Substituting the latter in Eq. J 1 8 b we get 

WU) = U ±±^f0L . (20) 

In order for Eq. J20b to be useful for our purposes, we shall apply 
some approximations, expressing it as a function of the phases 8. 
Appendix [C] also shows that for not unreasonably large values of 
the uncertainties, the first derivative of the perturbation function is 
f'(t e ) 1. Thus, it can be ignored in the denominator of Eq. i20i . 
The argument of the sine in Eq. < !20t has also to be expressed in 
terms of 9 = ut ew . Substituting Eq. J 1 5 b we obtain: 

9 = U)t ew = U)t e - w/(t e ) => t e = 9/0J + f(t e ), (21) 

and at a first order approximation we can set 

t e = 9/uj + f(8/cu). (22) 

Even if we do not give now a direct estimation of this approxima- 
tion, we shall realize that it is effectively safe when in Section 6 our 
analytical results will be compared with simulations. Putting it all 
together, the approximated pdf of the phases (9) is 

Pe(9) ~ U [l + sin(^0 + wo/(0/&;))] , (23) 

and substituting in Eq. d!4t t e by 9/u>, the perturbation function 
expressed in terms of Q/u> is equal to 

f(6/u) ~ 8k ■ sin ( ^-9 + <t> + V> J - 5c. (24) 



5 HOW THE PERTURBATION FUNCTION AFFECTS 
THE POWER SPECTRUM 

Once the pdf Pg of the phases perturbed by the errors on the orbital 
parameters has been evaluated, we can study how the power spec- 
trum is affected. In the Appendix lAlwe briefly introduce a method 
by which the power spectrum is directly expressed as a function of 
the pdf Pg. There, we find that a key role is played by the terms 

2V-1 

J2 Pb(2tti + 8), (25) 

i=0 
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which correspond to the pdf Pg folded in 2-k. Indeed, in Eq. J25b . 
^ 9 < 2k, the folding is due to the sum over the term 2iri, and 
N is the number of rotations made by the plausible neutron star 
during the whole observation T b s 



TV 



2-k 



(26) 



Figure [2] shows the folded Pg due to a 100% sinusoidal signal 
emitted by the pulsar, assuming two different sets of orbital param- 
eters. In the left panel, the correct set is used to calculate the Romer 
delay, so the perturbed function is null and 



\-Y^Pe{^i + 6) = l + sin0, 



N 



(27) 



where the factor 1/N cancels out setting U = 1 in Eq. d23 b . In 
the right panel, a value of the projected semi-major axis A slightly 
different from the true one is assumed. We can note that the folded 
Pg is still sinusoidal, but the effect of the wrong value of A is to 
reduce its amplitude so that 



N-l 

— J2 p e( 27li + 0) = 1 + e • sin(0 + , 



(28) 



Similarly, the power spectrum calculated at the signal frequency 
too is maximum when the right value of A is used (so when e = 1 
in Eq. l!28t). Using a different value of A, the power is reduced 
by a factor equal to e 2 (the square happens because of Eq. l |A5t 
in Appendix [Aj. This effect is not only related to the semi-major 
axis, but it is valid for all the orbital parameters, as we demonstrate 
below. 

At the signal frequency ujq, the term in the sum of Eq. l |28t is 
explicitly equal to 



Pg(2Ki + 9) = 1 + sin [9 + u ■ f 



2Ki + 9 



(29) 



This is obtained substituting the variable 9 with 2Ki + 9, and setting 
lj = wo in Eq. J23l >. The dependence of the perturbation function 
from 9 is negligible in Eq. l |29l >. Indeed, making the same substitu- 
tions in Eq. J24b we have 



'2ki + 9 
f ] — ok ■ sin 

w 



(30) 



Since here ^ 9 < 2-k, f2 or b/wo -C 1, and (r2 or b/cJo)# is al- 
ways very small, while the term (f2 or b/wo)27ri can be as large as 
^orbTobs = 27rT b s /P rb for i = N. Similarly, the dependence 
from 9 is negligible in Sk and ip. In what follows, the perturbation 
function will be labelled as fi, indicating that it depends only on 
the index i. 

With all this, the folded pdf becomes 



JV-l N— 1 



N 



(31) 



where the sum in the right hand is of a set of several sines with 
the same periodicity. The trigonometric identity of Eq. l[6j is a par- 
ticular case of a more general theorem stating that a sum of sines 
and cosines with equal periodicity, but different amplitudes and off- 
set phases is equal to a single sine with same periodicity, and with 
amplitude and offset phase depending on those in the sum. This im- 
plies that the sum in Eq. J3U is equal to e ■ sin(# + a), proving the 
equivalence with Eq. l |28l l. 

Summarizing: the effects of the errors of the orbital parame- 
ters (dp) on the power spectrum calculated at the signal frequency 



wo are described by the single factor e, so that P(wo,dp) = 
e 2 P(uj ,0). 

With a few steps of extra algebra we shall find a useful formula 
to evaluate the factor e 2 . The right hand side of Eq. < |28b can be 
written as 



1 + e ■ sin(# + a) — 1 + e [sin# cosa + cos9 sina] 
Similarly, the right hand side of Eq. d31t is equal to 

1 iV_1 

i=0 

1 iv_1 

1 + — Y, [sinS cos(wo/<) + cos6'sin(ajo/ l )] = 



(32) 



(33) 



i=0 
JV-l 



1 + sin^ — ^ cos(wo/i) + cos9 — sin(wo/^ 



N 4-i N 

!=0 

Comparing Eqs. l [32l and < f33b . we get 



cos(wo/i) = A^ecosa, 

i=0 
JV-l 

sin(ajo/i) = A^esina. 



And squaring and adding Eqs. J34b and l !35t we obtain 



Y sin(w /i) 



+ 



Y cos(w /i) 



N 2 e 2 . 



(34) 
(35) 

(36) 



6 CONSTRAINTS ON THE PARAMETERS 

So far we have commented on the way in which the perturbation 
function affects the power spectrum. In this Section, we aim to 
constrain the uncertainties (dp) in the parameters (p) in order to 
maintain the ability to detect pulsations. In practice, this reduces in 
searching for a formula that allows to state that if the uncertainty is 
smaller than a given value, dp < x, then e 2 > y. The larger is e 
the better, until for e = 1 there is no loss introduced by imprecise 
knowledge of the orbit. If one aims to search for pulsations from a 
compact object in a binary system, for which orbital parameters are 
known just to an indicative level, this study will provide the maxi- 
mum steps in the sampling so that the signal detection is secure at 
a certain level. 

Developing the squares of the two sums in Eq. i36l , it becomes 

^sin(w /i)] + [^cos(wo/i)] = 
N + 2cos(w /i)cos(wo/2) + 2sin(w /i)sin(w /2) + — = 

N + 2 [cos(w (/i - /a)) + ■•■] = (37) 



N + 2 



Y c ° s ( w o(/i - fi)) 



_i=l j — 



To continue further, we need now to discern the cases for which the 
partial derivative of the eccentric anomaly (dE/dp) is null, from 
those for which it is still a function of E. Indeed, in the former 
case (for the orbit parameters A, and W), Sk and ip in Eq. d30l > are 
constants that do not depend on the sub-index i of the perturbation 
function. 
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Pulse Phase (6) Pulse Phase (6) 

Figure 2. The histograms in both panels are the normalized folded phase distribution of a periodic signal emitted by a pulsar in a binary system with the 
parameters listed in Table [2] The emission is simulated as described in section 6.1.1, assuming a pure sinusoidal signal with a rate of 300 counts/day. The 
histograms are fitted by a sine function (black lines). Left: the true values of the orbital parameters are used to calculate the Romer delay. The amplitude of the 
sine function is equal to 1.0. Right: a value of the projection of the semi-major axis A different from the true one by 0.075 lt-s has been used to calculate the 
Romer delay. This causes the reduction of the amplitude of the sine function, that in this case is equal to 0.6. 



6.1 Cases for which 5k, and ip are constants 



imated as 



We consider the double sum in the square brackets of the last equa- 
tion. The difference of the perturbation functions (/» — fj) is equal 
to 



N-i 



. "u/^orb 



W /f2 o 



rb 



/] COS (w (/i - fi+An)) 



fi - fj = 



Wo 



N-i 



Wo /fiorb 27rf2 or b J 



i (too ■ 5k[Ki — sin(Ai + x)]) dx 



= 8k 



I — — 27ri + <f> + tp] — sin ( — — 2irj + cj> + tp 



Wo 



Wo 



N-i 
2tt 



Ct, (41) 



5k 



Ki — sin 
= 8k 



W 



2n(i + An) + <jj + tp 



(38) 



Ki — sin ( ^° rb 2n An + Aj 

U) 



where Ki has values between —1 and 1, and we set j — i + An. In 
this way, for each fixed value of the index i, the difference (/; — fj ) 
is function of An, and the sum over j in Eq. J37t becomes 



^ COS (wo(/i - fi+An)) ■ 



(39) 



Furthermore, (fi — fi+An) is periodic, with period An = 
ojo/f2 rb- This feature allows us to solve, with a good approxima- 
tion, Eq. d36t . and evaluate e 2 . 

The argument of the sum in Eq. l |39t has the same periodicity 
of (fi — fi+An)- Then, we can approximate it as 



/] COS (wo(/< - fi+An)) 



^ COS (wo(/t - fi+An)) ■ (40) 



N-i 

W0 / fiorl 



In the right hand side of Eq. J40b the sum is over one cycle, while 
the term (N — i)/(wo/f2 rb) is the number of cycles of the func- 
tion (fi — fi+An)- This approximation is good when the number of 
cycles is large. Actually, we will show that it is still good when the 
full observation includes just one cycle or more, whereas it starts to 
be annoyingly imprecise when less than one cycle is observed. 

Setting x = An ■ 27rfi or b/wo, Eq. J40t can be further approx- 



where the integral has been referred to as d, indicating that it de- 
pends only on the index i. Note that since the integral is over one 
cycle (0, 27r), the phase Aj does not have any effect, so we can 
safely set Aj = 0. Then, d is a function of the term Ki, that we 
already noticed has —1 and 1 as its minimum and maximum values, 
respectively. Considering Ki as a continuous function, the average 
value of d is equal to the double integral 

d = - / cos (wo ■ 8k[Ki — sin a;]) dxdKi. (42) 
2 J-i Jo 

The order of the integrals can be inverted, and an analytical solution 
for the integral on dA', can be found. 

j cos (ujoSk ■ [Ki — sin a;]) dKi 
sin(wo5fc • Ki) ■ cos(wo<5fc • sin a;) 



ujoSk 



(43) 



(44) 



=> j cos (wo<5fc • [Ki — sinx]) dKi = 

2sin(wo5fc) • cos(wo<5fc ■ sin a;) 
wo<5fc 

Finally, the average value of C\ is equal to 



— sin(w 5fc) f 2w , . . , .... 

d = — / cos(woo/c • smx)dx. (45) 

woofc Jo 

Summarizing, the sum in Eq. i39\ is equal to Eq. d4 H : 

N-i 



^ COS (wo(/» — fi+An)) = ~Ci. 



2tt 



(46) 
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Table 2. Parameters of the simulated pulsar. 



"_ -— eq.36, T oba = P OIb ; T obs = 10p o ,b 
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— eq.50 / M \ \ 









-0.5 0.5 

<o n 5k/ji 



Figure 3. The factor e 2 calculated with different approximations in the case 
described in Section 6.1. See the text for the details. 



Substituting it in Eq. 1371 l. and then in Eq. d36t we get 



JVV = N + 2 



,i=l 



2tt 



(47) 



Taking into account the average value of d, this equation approx- 
imately becomes 



N 2 e 2 



N + 2 



= N + 2 



' N 

.t=i 

N 2 - N 



2tt 



271 



2tt V 2tt 



(48) 



In conclusion, we have found that for N ^> 1, the factor e that 
measures the loss in the power spectrum when values of the orbital 
parameters are not precisely known is approximately 

2 1 sin(wo(5fc) f 2n , „ 

e = — ■ — — / cosfojoofc ■ sm x)dx. (49) 

27T uj 5k J 

This result gives a clear description of the power spectrum 
around 5 k — 0, from which the following conclusions can be de- 
rived: 

• The factor e 2 is function only of ujoSk. The other parameters 
(4>, ip, 5c) do not affect the power spectrum. 

• Eq. d49t does not depend on the duration of the observa- 
tion T bs, at least within our approximations. Figure [3] shows 
the factor e 2 as a function of ujoSk, as expressed in Eq. ( I49t - 
Also in Figure [3] we show the values of e 2 calculated by the 
non-approximated formula of Eq. i36\ . for T b s = 10P or b, 
Tobs = Parb, Tobs = 0.5P or b- As we can notice from the 
plot, for T obs > P orb , Eq. (@9j is a good approximation of e , 
whereas it underestimates the values of e 2 for T b s < Porb- This 
feature is particularly advantageous for pulsation searches. Indeed, 
the maximum errors on the orbital parameters needed to avoid 
washing out the periodic signal are unaffected by the duration of 
the observation, which, on the other hand, the longer it is, the 
higher is the signal to noise ratio of the power spectrum, which 
is proportional to the total number of events (Scargle 1982). We 
can deduce from Figure [3] that reducing the observation time to 
Pobs <S Porb, the factor e 2 remains unaffected. This is what most 
commonly happen in radio observations, which are so short that 



Parameters 


value 


p 

1 psr 


300 ms 


Porb 


4.0 days 


T 


54587.00 MJD 


A 


2.5 lt-s 


W 


44.39° 


e 


0.61 


T^obs 


40 days 



the orbital motion has negligible effects on the pulsation search. 

• The integral in Eq. l !49t can not be solved analytically, but 
its behavior is very similar to the sine function. Indeed, we have 
empirically found that a very good approximation for e 2 is 



(50) 



2 

e = 




2 













As shown in Figure [3] Eq. d50t is even better than Eq. l |49t to de- 
scribe the central peak of e 2 , but it fails in the side lobes. 

Clearly, these considerations are valid for the case analyzed 
in this section, i.e., when the parameters 5k and ip have no depen- 
dence on the eccentric anomaly E. This happens when the deriva- 
tive of the eccentric anomaly with respect to the canonical orbital 
parameter (dE/dp) is null, as is the case for the projection of the 
semi-major axis A, and the longitude of periastron W. 



6.1.1 Semi-major axis A and simulations 

When we take into account the error on A, uioSk is explicitly equal 
to(seeEqs.[8]-[LQ 



uj Sk = uj dA^/l-e 2 cos(W) 2 . 



(51) 



In order to maintain the factor e 2 higher than a certain level (say 
e 2 > 0.4) the term ujoSk/TY has to be lower than the inverse value 
(u) 8k/iv < [e 2 ] _1 (0.4) ~ 0.4), meaning that 



dA < 



oj ^l-e 2 cos(VT) 2 



(52) 



Clearly, the error dA is inversely proportional to the pulsar fre- 
quency (wn). It is interesting to notice that high eccentricities are 
less constraining for dA, even though the term \Jl — e 2 cos(W) 2 
change very slowly with the eccentricity. For example it is lower 
than 0.5 only for e > 0.85, and cos(W) ~ 1. 

In order to check these results, we simulated the barycentred 
arrival time series from a pulsar in a binary system with the features 
described in Table|2] Then we demodulate it, but modifying one of 
the orbital parameters (in this case A). In this way we obtain the 
demodulated arrival time series perturbed by the variation of the 
orbital parameter, called hereafter perturbed time series. Its power 
spectrum is expected to follow the analytical description described 
above. 

The simulation consists of two steps. First, the emission time 
series in the pulsar reference frame has been created. All the time 
stamps are taken as random numbers in the range [0; 27r/wo], fol- 
lowing a pure sinusoidal distribution 



Pt{t) = l + sin(o;o<), 



(53) 
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where ujo is the assumed frequency of the pulsar. Then, in order to 
cover the full duration of the observation, each time stamp has been 
randomly delayed adding a value {2-K/uJo)n, where n is a random 
integer number uniformly distributed in [0, woTobs/27r]. Finally, to 
pass from the emission times in the pulsar reference frame (t e ), to 
the barycentred arrival time (t a ), the Romer delay An(t e ) is added 
to each t e . The Romer delay is calculated taking into account the 
orbital parameters. The demodulation is p erformed by means of the 
algorithm describe d inlRea et alj d2011a|). which makes us e of the 
program TEMP02 l lHobbs. Edwards. & Manchester! ( 120061) ). 

The power spectra shown in Figure|4]are calculated from time 
series demodulated when varying the projection of the semi-major 
axis by three different amounts dA. In the top panel dA = 0.053 
lt-s correspond to cooSk/n — 1/n. From Figure[3]or Eq. \50\ , for 
this value of LOoSk/ir the power calculated at the pulsar frequency 
loo is expected to be ~ 0.6 times the coherent power resulting from 
an unperturbed demodulation. This is confirmed by the power spec- 
trum in the top panel in figure|4] In the middle panel dA — 0.067 
lt-s, corresponding to ujoSk/n = 0.4, and e 2 ~ 0.4. In the bot- 
tom panel dA — 0.125 lt-s, implying uoSk/n — 0.75, for which a 
null power is expected. Indeed, the peak at the pulsar frequency is 
very suppressed, while sidebands dominate in this case. The ratio 
of the power to the coherent power of the central peak in the plots 
of Figure|4]can be directly compared with the e 2 evaluated in Fig- 
ure[3]at the corresponding values of cuoSk/n. We thus find a good 
agreement between the analytical procedure and the simulations. 



6.1.2 Sidebands ambiguity and a method for removing it 

As shown in the bottom panel of Figure [4] it is possible that for a 
certain value of dA, the peak at the true pulsar frequency is strongly 
suppressed, while the highest sideband peak is detected above the 
noise level. How can we realize in the course of a real observation 
whether the peak detected is a sideband or the pulsar itself? Be- 
fore answering to this question, it is important to clarify that the 
presence of the sideband structure in the power spectrum is a clear 
signature of the presence of a pulsar in the system, since otherwise 
we would see only n oise. 

iRansoml l l200ll) studied the sidebands by a not-demodulated 
arrival time series. The number and the height of the sideband peaks 
depend on the orbital parameters. Here we summarize some impor- 
tant features of the sidebands, which can help solving the ambiguity 
described above. 

• The difference in frequency of two adjacent peaks is equal to 
the orbital frequency, so that with respect to the pulsar frequency, 
the sideband peaks are located at 



ojsb = Wo ± ntt 



(54) 



where n is an integer number, and uj s b the frequency of the n'th 
sideband peak. 

• The sum of the power of all the sideband peaks is equal to the 
total coherent power of the pulsar. 

• The larger is the orbit, more numerous are the sideband peaks, 
because the modulation in the arrival time series is stronger. 

In our case, the sidebands appear in the power spectra of the 
perturbed time series because they still have a residual modulation. 
This means that larger is dA, the higher the residual modulation 
will be, and because of the third bullet commented above, more nu- 
merous sideband peaks will appear. We can deduce that the highest 



probability to detect a sideband peak, rather than the pulsar one 
is when e 2 is close to the minimum in Figure [3] Indeed, when it 
happens, the pulsar peak is almost at zero power and consequently 
the sideband peaks get stronger. In this case, their number is small 
because we are close to the true value of the orbital parameters. 
Therefore, at least one sideband peak could have a power higher 
than the noise level. 

Thus one can think of a method to solve the ambiguity intro- 
duced by the possible presence of sidebands as follows; 

(i) When a significant peak is detected at in the power 
spectrum, we can assume it to be at the minimum of the e 2 curve, 
regardless on whether this is true or not. 

(ii) With this assumption we can get an estimation of dA (or 
more in general of Sk) by means of Eq. l !52t setting in it e 2 = 0, 
and cjo = uid (since Q OI b <C wo, for Eq. l |54t it is always true that 

U) d — Wo). 

(iii) We can use the so-estimated dA to define a sub-sampling of 
the parameter space in order to explore the profile in Figure [3] and 
to understand if we really are in its minimum, or on the top of the 
peak. Using the new sampling, and if in the former case, it is ex- 
pected that a frequency peak with higher power with respect to the 
detected one will appear at a frequency oj = aid ± nfiorb- This new 
peak is then the true pulsar frequency. In contrast, if the new peak 
does not appear, or its power is lower than the first detected one, 
the true pulsar frequency remains wq = tod, as originally detected. 

Since the presence of the sidebands in the power spectrum 
means that there is a pulsar in the system, their detection will imme- 
diately solve the nature of the compact object if such is unknown. 
A techni que dedi c ated t o the dete ction of sidebands w as formu- 
lated by IRansoml d200ll) (see also iRansom et all l2003h . In brief, 
because of the first feature listed above, the sideband peaks appear 
as a short series of regular pulsations in the power spectrum of the 
arrival times series. Then, the detection of the sidebands is possible 
by taking the Fourier transform of this short section of the power 
spectrum, which is expected to have a peak at the orbital frequency 
of the system. This technique was originally formulated for the de- 
tection of radio binary systems with short period P or b < T b s - 
Since this condition is common in gamma-ray observations, this 
technique can be easily adapted to the perturbed time series studied 
in this work. 

Whether the first or the second method described above is 
more appropriate, strongly depends on the signal to noise ratio, as 
well as on the pulsed fraction of the signal. For weak pulsed signals 
the first method should be more appropriate, because most of the 
sideband peaks would be lower than the noise level. On the other 
hand, when the method proposed bv lRansonj fcOOlh is applied to 
calculate the significance of the signal one should take into account 
also the trials in Fourier transforming several short section of the 
power spectrum. 



6.1.3 Longitude of the periastron W and simulations 

For the longitude of the periastron, the expression of uioSk is 



r /l- e 2 + e 4 C os 2 Wsin 2 VK 
D 8k = uj AdW\ — — . (55) 

1 — e 2 cos 2 W 

The error dW is inversely proportional to both the pulsar frequency 
and the projection of the semi-major axis A. Its dependence on W 
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Figure 4. Power spectra of the simulated time series, demodulated varying the projection of the semi-major axis by three different amounts dA. Top panel: 
dA = 0.053 lt-s correspond to uj()Sk/n = 1/tt, and e 2 ~ 0.6 is expected (dashed line). Middle panel: dA = 0.067 lt-s, correspond to uioSk/ir = 0.4, and 
to an expected e 2 ~ 0.4 (dashed line). Bottom panel: dA = 0.125 lt-s , imply iooSk/n = 0.75, a null power is expected. The so-called sidebands dominate 
in this case. The Coherent Power at the denominator of the y axis is the power calculated for an unperturbed demodulation (dA = 0) at the pulsar frequency 

U1Q. 



and e is in the square root. For W = n, this term is equal to 1, 
so it does not affect dW. In contrast, for W = tt/2, or W — 
|-7r the square root is equal to yl — e 2 , so also in this case high 
eccentricities are less constraining for dW. 



We have demodulated the arrival time series of the pulsar 
simulated in the previous Section by varying W in three different 
amounts, dW. Figure [5] shows the power spectra of the perturbed 
time series. In the top panel dW = 1.34° implies un)8k = 1, 
and e 2 ~ 0.6 by Eq. [50] In the middle panel dW = 1.68°, 
cooSk/iv = 0.4, and e 2 ~ 0.4 (Eq. [50] see also Figure 0. In the 
bottom panel dW = 3.15°, uj 8k/n = 0.75, and e 2 ~ 0. Also 
in this case the sidebands dominate in the bottom panel, while the 
peak at the pulsar frequency is prominent in the top one. 



6.2 Cases for which Sk and ip are functions of the eccentric 
anomaly E 

Similarly to the previous Section, we would like to find a solution 
for e 2 when the parameter known with the higher uncertainty is 
the eccentricity e, or the epoch of the periastron To, or the orbital 
period P or b- These cases are more complicated because the partial 
derivatives of the eccentric anomaly E with respect to these pa- 
rameters are not null. This implies that Sk and ip are not constants, 
but rather, functions of the eccentric anomaly E. In the following 
subsections we analyze these three cases one by one. 



6.2.1 Eccentricity e and simulations 

The differential of Eq. d!3t with respect to the eccentricity e leads 
to the following expression for the partial derivative of E respect 
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Figure 5. Power spectra of the simulated time series, demodulated varying the longitude of the periastron by three different amounts dW. Top panel: dW = 
1.34° correspond to u)Q8k/n = 1/tt, and e 2 ~ 0.6 is expected (dashed line). Middle panel: dW = 1.68°, correspond to uJoSk/n = 0.4, and to an expected 
e 2 ~ 0.4 (dashed line). Bottom panel: dW = 3.15°, imply uioSk/n = 0.75, a null power is expected. The Coherent Power at the denominator of the y axis 
is the power calculated for an unperturbed demodulation (dW = 0) at the pulsar frequency ujq . 

to the eccentricity 

sin(f2 or bt e ) 



Note that Ski and ipi contain the partial derivative of E as function 
of (2m + 9) /ojo, which is equal to 



dE 



sinE 



de 1 — ecosS 1 — e cos($l or bi e ) 



(56) 



Since dE/de is periodic, so are Sk and ip. In order to find a solution 
for s 2 , we need to solve the double sum in the square brackets of 
Eq. l !37t . where fi — fj is in this case 



dE /2ni + i 
de \ oj 



■ "o "n ' 

1 - e cos (Sari. 27ri + Sari 

%^2tti) 



■ ecos(^f-27ri) 



(58) 



fi - fj = Ski ■ sin 



a, 



wo 



-2m + <(> + - (57) 



Ski • sin ( — ^-2nj + (j> + ipj 



wo 



Ki — 8ki + A n ■ sin 



UJ 



■2ty(i + An) + (/> + Tpi+An 



Since this is periodic in Q OI h/uo, then also f% — fj = f% — fi+An 
is periodic with period An = fl O zb/^0- This is the same fea- 
ture we have found in the previous Section. Hence, we can follow 
the same reasoning and integrate the argument of the double sum 
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(cos(uuo(fi — fj))) in one cycle, as we did in Eq. J4U . to get 

N-i 



E 



)S (w (/i - fi+An)) 



N-i 



u 



f 

Jo 



cos (uio ■ [Ki 



wo/n or b 27r51 or b 
8k(x) ■ sin(x + ip(x))]) dx 
N 



2tt 



-a 



(59) 



Here _K"i = Ski ■ sin ^ 2-7ri + + ipi^J is a periodic function 
that oscillates between the values K m i n and K max . Following the 
reasoning of the previous Section, from Eqs. d42t and < l49t > we have 
in this case that 



2 

e = 



27r(A"„ 



K„ 



x 



cos (wo ■ [Ki — Sk(x) ■ sin(x + ip(x))]) dxdKi 



(60) 



2 

e = 



The integral over dKi can be analytically solved, obtaining 

1 sm(u!oK max ) — sm(ujoK 



■2tt 



uj (K v 



-K r 



cos (ojo ■ Sk(x) ■ sin(x + ip(x))]) dx. 



(61) 



Here the values K m in and K max are not symmetric. Hence the ex- 
pression outside the integral can not be further simplified to a sine 
function. Furthermore, Kmin and K max cannot be expressed by an 
analytic formula. However, some features can be found analyzing 
the function Ki : 



Ki = Ski • sin 



a-, 



+ 



u 
Ade ( ( 
etan(W) 



2m + <j> + ipij = 
e 2 cos 2 W)x 



■ e 2 + tan(W) ^T- 
e 2 cos i W 



+ 



1/2 



sin x + atan 



/ tan(W) 



+ atan — 



e 2 cos 2 W / 
1 - e 2 cos 2 W 



(62) 



1 



+ 



etan(W) 



e cos 2 W 
1 



l-e 2 +tan(W0 ^1^2 



where we have set x — (f^ or b/a;o)27ri. Eq 
read as Ki = A deg(x; e, W), so that K 



can be more easily 

= Adeg max and 



K m i n = Ade g-min ■ Substituting in Eq. ( 1611 ) it is clear that ojo A de 
has a key role in this case, from which we can conclude that the 
maximum error of the eccentricity (de) that is allowed to avoid 
washing out the pulsed signal is inversely proportional to the pulsar 
frequency (cjo), and to the projection of the semi-major axis (A). 

Also in this case, Eq. i ]6 1 | i does not depend on the duration 
of the observation (T b s ), when it is greater than one orbital pe- 
riod (P or b). In order to check this, Figure [6] shows Eq. l |61| > versus 
coo Ade for different values of e and W. They are compared with 
Eq. \36l assuming different observation times (T b s = 10P O rb, 
Porb, and 0.1P or b). From these plots we can notice that the agree- 
ment between Eq. d36t and Eq. d6U is not always good. Further- 
more, the width of central the peak produced by Eq. J36b , and its 
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Figure 6. e 2 versus ui^Adej-K as calculated using Eq. {63, and Eq. (36) 
for the case described in Section 6.2. 1 . Each panel assumes different values 
of e and W. Top: e = 0.3, W = it. Middle: e = 0.7, W = w. Bottom: 
e = 0.7, W = 1.47T. Three different observation times are assumed to 
compute Eq. (36): T b s = 10P or b; -P rb. al 'd 0.1P or b- 



first zeroes have an evident dependence on e and W. We investi- 
gated this dependence assuming a grid of values of e and W, and 
numerically finding the values of ujq A de/ir for which e 2 = 0.4. 
They are represented by the surface plot of Figure [7] (top panel), 
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Figure 7. Top: Surface representing the values of u)oAde/iT for which 
e 2 = 0.4 versus e and W. This surface is approximately fitted by the 
function S(e, W) given in Eq. i63i . Bottom: The almost-flat surface in this 
plot represents the values of ljq Ade/nS(e, W) for which e 2 = 0.4. The 
flatness of the surface indicate that the term LUoAde / 7rS(e, W) absorbs the 
dependence of e 2 on e and W . 



which is approximately fitted by the function 



S(e, W) = (7(e) - Ampl(e) ■ \cosW\ , 

U(e) = U + a • Log(6 - 
Ampl(e) = Ampl • (1 - e/2)e 2 



with 



(63) 



a 1 ' 3 ), 



and where U = 0.816 ± 0.001, a = 0.345 ± 0.005, b = 1.08 ± 
0.01, and Amplo = 0.84±0.01. If we now plot the previous values 
of LJoAde divided by S(e, W), we obtain an almost flat surface at 
z = 1 with small fluctuations due to the approximations (see Figure 
|7] bottom panel). The flatness of the surface indicate that the term 
uj()Ade/TvS(e, W) absorbs the dependence of e 2 on e and W . In 
other words, the factor e 2 is approximately a function of the single 
variable z = ujq Ade/nS. To show this better, Figure [8] shows the 
values of e 2 obtained using Eq. l |36t versus DoAde/nS, assuming 
for e and W the same values as for the three plots in Figure [6] 

We can now state that in order to maintain the factor e 2 higher 
than a certain level, the following condition must be satisfied 



ojq Ade 



< [e 2 ]' 1 ■ Me) + Ampl(e) ■ jcosW|} . 



(64) 



Figure [9] shows the functions 17(e) and Ampl(e) fitted on a set of 
values of e. Looking at this figure, and from Eq. l !64t . we can con- 
clude that for low eccentricities (e < 0.3) Ampl(e) is negligible, 



Figure 8. e 2 versus uoAde/irS, assuming for e and W the same values as 
for the three plots in Figure[6] The values of e 2 are calculated by means of 
Eq. {36). Blue line: e = 0.3, W = n. Red line: e = 0.7, W = n. Green 
line: e = 0.7, W = 1.4tt. 



and U (e) acts like a weak scale factor ~ 0.8. In contrast, for high 
eccentricities (e > 0.6) and especially for W ~ 1, the error de is 
strongly constrained if we wish to maintain pulse-detection capa- 
bility. Indeed, the right hand of Eq. J64t can reach values lower than 
0.2 for e > 0.9, implying that de has to be lower than 0.2-7r/(u;oA) 
in order to avoid a prohibitive suppression of the pulsed signal 
(e 2 < 0.4). 

To check on all these results we used the arrival time series 
of the pulsar simulated in Section 6. 1 . 1 . The demodulation is per- 
formed varying the eccentricity by three different amounts de, and 
the corresponding power spectra are shown in FigureQJJ In the top 
panel de = 0.0107 correspond to ujuAde/nS = 1/n. In the mid- 
dle panel de = 0.0337, corresponding to ujoAde/nS = 1, and 
e 2 = 0.4 as discussed above. In the bottom panel de = 0.0675, 
implying uioAde/nS = 2, for which a null power is expected as 
we can deduce from Figure[8] The ratio of the power to the coherent 
power of the central peak in the plots of Figure[l0]can again be di- 
rectly compared with the e 2 as computationally evaluated in Figure 
[8]at the corresponding values of uioAde/nS . The very good agree- 
ment between the analytical and the simulation results validate the 
procedure adopted for the eccentricity. 



6.2.2 Epoch of the periastron To and simulations 

The treatment of the uncertainty of the epoch of the periastron is 
similar to the previous case. From Eq. d 1 3b . the partial derivative of 
the eccentric anomaly E with respect to To is 



dE 



rb 



e cosT 



(65) 



This is periodic, and implies that also 5k and ip are. In particu- 
lar, the terms Ski and ipi in Eq. l !57t contain the partial derivative 
of E as function of (27ri + 6)/ ojq, which is periodic with period 
fiorb/wo- Since this is the same feature we have found in the case 
of the eccentricity, we can follow all the steps done in the previous 
Section until Eq. <6U . 

Also in this case the values K m i n and K max are not sym- 
metric. Hence the expression outside the integral can not be further 
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Figure 10. Power spectra of the simulated time series, demodulated varying the eccentricity by three different amounts de. Top panel: de = 0.0107 corre- 
sponds to uioAde/nS = The horizontal dashed line indicate the expected value of e . Middle panel: de = 0.0337, corresponds to uioAde/nS = 1, and 
e 2 = 0.4 is expected (dashed line). Bottom panel: de = 0.0675, implies uJoAde/nS = 2, a null power is expected. The Coherent Power at the denominator 
of the y axis is the power calculated for an unperturbed demodulation (de = 0) at the pulsar frequency uiq. 



simplified to a sine function. However, the function Ki is equal to: 



Ki — Ski ■ sin ( ^ orb 2ivi + <p + ipi 



fiorbAVl - e 2 cos 2 W 



dT o -sin(s: + 0-7r/2) 



QovhA 

1 + e 2 



■ (1 + e Vl ~ e 2 cos 2 WdT • £2!^± 



1 — e cosx 
adT ■ g(x; e, W) (66) 



where we have set x — (f2 or b /uJo)2iri, a = O or b A/ (I + e 2 ), and 
g{x;e, W) = (1 + e 2 )Vl — e 2 cos 2 W cos(a; + cj>)/(l — e cosa;). 
We have found useful to introduce the term (1 — e 2 ) in Eq. ( |66| > for 
the comparison of the analytical results with several of the simula- 
tions we performed. Substituting K max and K m i n in Eq. d61t with 
the corresponding maximum and minimum of Eq. j66\ , it is clear 
that LuoadTo has a key role in this case, from which we can con- 
clude that the maximum error of the epoch of the periastron (To) 



that is allowed to avoid washing out the pulsed signal is inversely 
proportional to the pulsar frequency (luq), and to the projection of 
the semi-major axis (A), while it is directly proportional to the or- 
bital period P orb = 27r/Sl orb . 



The panels in Figure QT] show e 2 versus ujoadTo for different 
values of e and W. The values of e 2 are calculated by Eq. ( |61t . and 
by Eq. ( 136b for different observation periods . The plots show that 
there is a good agreement between Eq. d6U and Eq. ( 1361 ). We note 
again that e 2 does not depend on the duration of the observation 
(Tobs) when it is greater than one orbital period (P or b); but the 
width of the central peak still has a dependence on e and W . To 
normalize it, we have used the same approach adopted in the case 
of the eccentricity. Namely, we have found the surface of uioadTo 
for which e 2 = 0.4 varying e and W, and we have searched for an 
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Figure 9. U(e) and Ampl(e) evaluated for several values of the eccentric- 
ity, and fitted by the second (red line) and third (blue line) functions of the 
group of Eqs. {63). 



approximate fitting function, that in this case is 
Ampl(e) 



U(e) 
Uo + a ■ Log(6 - 
Ampl • (1 - - 



(l + cos(2W0), with(67) 



:) e 



S(e,W) 

U(e) 
Ampl(e) 



where U = 0.410 ± 0.003, a = 0.25 ± 0.01, b = 1.07 ± 0.01, 
Amplo = 0.45 ± 0.01, and c = 0.68 ± 0.01. In this way we have 
found that the factor e 2 is approximately a function of the single 
variable z = U)oatdTo/irS. This is evident in Figure[T2] where e 2 
obtained using Eq. j36\ is plotted versus uootdTo/nS, assuming 
for e and W the same values as in Figure [TTI 

In FiguresQT](bottom panel) and[l2](green curve) we can no- 
tice that for W ~ n the curves of e 2 obtained by Eq. ([36} do not 
reach zero in the neighbourhood of ujoadTo / n S = 0. In particu- 
lar, when e > 0.9, the first local minimum closer to the peak of the 
curve has a value e 2 > 0.4. From the simulations we performed 
we observed that this feature is true also for values ofWj^ir. Un- 
fortunately, this feature is lost in the approximations done to pass 
from Eq. l |36t to Eq. \6l\ , as shown in Figure [TT1 (black curves). 
Neglecting this, we can conclude that in order to maintain the fac- 
tor e 2 higher than a certain level, the following condition must be 
satisfied 

WoctdTo 



< £ 



U(e) 



Ampl(e) 



(l + cos(2W)U, (68) 



where the inverse function [e 2 ] -1 can be deduced by Figure [72l 

To check on these results we have also used the arrival time 
series of the pulsar simulated in Section 6. 1 . 1 . The demodulation 
is performed varying the epoch of the periastron by three different 
amounts dTo, and the corresponding power spectra are shown in 
Figure [T3l In the top panel dTo = 5.2 x 10~ 3 days corresponds to 
uoadTo /ttS = 1 /it. In the middle panel dT = 1.63 x 10~ 2 days, 
corresponding to ujoadTo/irS = 1, and e 2 = 0.4 as discussed 
above. In the bottom panel dTo = 3.26 x 10~ 2 days, implying 
LJoadTo/irS = 2, for which a null power is expected. The ratio 
of the power to the coherent power of the central peak in the plots 
of Figure [T3] is in good agreement with the values of e 2 as com- 
putationally evaluated in Figure [12] at the corresponding values of 
uoadTo / it S , with the exception of the bottom panel in Figure [T3l 
where the central peak has a power of ~ 0.2 against the null value 




eq.36, T obs = P ort ; 10 

et "- 36 ' T .b,= "- 1P .rb 

eq.61 



0.7; W = 71 




-2 -1.5 



Figure 11. e 2 versus LOoadTo /ir as calculated using Eq. {6TJ. and Eq. {36) 
for the case described in Section 6.2.2. Each panel assumes different values 
of e and W. Top: e = 0.3, W = 0.6tt. Middle: e = 0.7, W = 0.6tt. Bot- 
tom: e = 0.7, W = 1.07T. Three different observation times are assumed 
to compute Eq. (36): T Q b s = 10P or b; -P rb> an d 0.1P or b. 



expected. This discrepancy is due to the neglected dependence of 
the level of the first local minimum of e 2 by the eccentricity, as 
explained above. 
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Figure 13. Power spectra of the simulated time series, demodulated varying the epoch of the periastron by three different amounts dTo . Top panel: dTo = 
5.2 X 10 -3 days corresponds to uqckITq j 'tt S = 1/tt. The horizontal dashed line indicate the expected value of e 2 . Middle panel: dTo = 1-63 X 10 — 2 
days, corresponds to uJoadTo /nS = 1, and e 2 = 0.4 is expected (dashed line). Bottom panel: dTo = 3.26 X 10 — 2 days, implies LuoadTo / nS = 2. The 
Coherent Power at the denominator of the y axis is the power calculated for an unperturbed demodulation (dTo = 0) at the pulsar frequency ujq ■ 



6.2.3 Solution for orbital period P or b and simulations 

From Eq. d!3l >. the partial derivative of the eccentric anomaly E 
with respect to P or b is 



BE 

dPorl 



dE dSl or b 



27r(l - ecosE) ~ 

Oorb-E 1 

' _ 2tt(1 - ecosE)' 



(69) 



Eq. l !69t is not periodic. Therefore, the procedure followed in the 
previous cases can not be applied for the orbital period. Eq. l |69t 
has instead a linear dependence on t e — To- As we will show, this 
implies that the effect of the error <fP or b does depend on the dura- 
tion of the observation T b s , in contrast with the other parameter 
errors. Substituting Eq. l !69t in Eqs. Q It . J 1 2b . and l |30t , we obtain 
the perturbation function due to the uncertainty of the orbital period 



fi = -A 



dPorb n o ^ttfi cosix + i 

— — vl — e 2 cos 2 VK 



(70) 



where x = ^jj^-2ni. The factor e 2 can be calculated substituting 
Eq. l !70t in Eq. J36t . Figure [T4l shows e 2 versus cjoAdP rb/Porb 
calculated for three different observation times (T b s = lPorb, 
T b s = 2P or b, and T b s = 10P or b), and random values of the 
e and W (e = 0.5, and W = 7r/3). The dependence of e 2 on T b s 
is evident in Figure Q5] We found that close to the peak, the factor 
e 2 is roughly approximated by 



i n 2 / a dPoib 
Wn I uiqA— 



(71) 



where n = Pobs/Porb is the number of orbital periods observed. 
Applying Eq. J7U to the pulsar simulated in Section 6.1.1 (see Ta- 
ble O, a decrease of the power by a factor e 2 = 0.9 is expected 
if the arrival time series are demodulated varying the orbital period 
by dP OI b = 7.6 x 10~ 4 days. This is confirmed in Figure [T5l 

In conclusion, the condition to maintain e 2 higher than a cer- 
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Figure 12. e 2 versus uocedTo / it S ', assuming for e and W the same values 
as for the three plots in Figure [T71 The values of e 2 are calculated by means 
of Eq. (36). Red line: e = 0.3, W = Q.6n. Blue line: e = 0.7, W = 0.6tt. 
Green line: e = 0.7, W = 7r. 
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Figure 14. e 2 versus uJoAdP or ^/P or ^, assuming e = 0.5 and W = 7r/3. 
The values of e 2 are calculated using Eq. j36\ . Thi'ee different observation 
times are assumed: T obs = P or b (blue line), T obs = 2P orb (green line), 
T bs = 10P or b (red Une). 

tain level (say e 2 > 0.9) can be derived by Eq. d71l > as 



Ujl orb ± I ± & 

Porb mj A\l 10 

Eq. d72t shows that the maximum error dP or b useful to maintain 
the loss in power e 2 above a certain value is inversely proportional 
to the pulsar frequency (wo), to the projection of the semi-major 
axis A, and also to the lenght of the observation T b s = nPorb- 



7 DISCUSSION AND CONCLUSIONS 

We have presented an analytical study aimed to understand the im- 
pact of the uncertainties of the orbital parameters on the pulsation 
searches. We validated the analytical study with numerical simula- 
tions. We especially focused on the cases where the observations 
span a time longer than the orbital period of the system. This is the 
usual case, for instance, in GeV observations with the Fermi-LAT 
gamma-ray satellite. 



The search for pulsations is performed by calculating the 
power spectrum of the emission times, which are evaluated demod- 
ulating the arrival time series. But the errors on the orbital param- 
eters unavoidably leads to a wrong estimation of the Romer de- 
lay, and consequently to an uncertain demodulation. Here, we have 
studied the relation between the exact and the imperfect emission 
time series considering that the latter are a perturbation of the for- 
mer. Within this frame, we have defined the perturbation function 
(Eq.ll4ll. and have found how it affects the probability density dis- 
tribution of the phases attributed to each event time assuming a 
generic frequency (Eq.l23ll. 

The power spectrum calculated at the pulsar frequency is re- 
duced by a factor e 2 when the demodulation is not correct. Eq. ( 1361 ) 
describes the relation among the perturbation function and the fac- 
tor e 2 . Starting from Eq. ( 136b . we have analyzed, in a case by case 
basis, the impact of the uncertainties of each of the orbital parame- 
ters. FigurefJ] and Eq. d50t show the behavior of the power loss, e 2 , 
in the cases of the semi-major axis A and the longitude of the peri- 
astron W. The maximum value of the error on the semi-major axis 
A in order to maintain the factor e 2 larger than a certain amount is 
explicitly given by Eq. d52| l. while in the case of W it can be de- 
duced from Eq. d55| l. The maximum error for the eccentricity and 
the epoch of the periastron are given by Eq. d64t and Eq. {68}, re- 
spectively. Whereas for the orbital period, the maximum error can 
be estimated by Eq. d72| >. All this is summarized in Table[3] 

The results discussed so far concern the case in which only a 
single parameter is affected by a significant uncertainty. The total 
loss factor that reduce the power of the coherent signal is expected 
to be the product of the factors calculated for each single parameter: 

2/ 2 2 2 2 2 

e — ea ■ £w • e e ■ Et ■ £ p orb ■ \i->) 

To understand this suppose that only the value of A differs from 
the true one. The measured power at the pulsar frequency is then 
reduced by P(cjo, dA) = £ A P(ljq, 0). If from this starting condi- 
tion the value of W is varied by dW, we expect a further reduc- 
tion of the power, such that P(ujo, dA, dW) = e 2 v P(a;o, dA) = 
e w ' SaP^o, 0). This reasoning leads to Eq. d73l l. Anyway, be- 
ing this only an intuitive demonstration, we can not exclude that 
combinations of the uncertainties on different parameters can devi- 
ate from the expected behavior of the total factor e 2 . In this sense, 
Eq. d73t . and the results that follow in this discussion have to be 
considered as conservative. In particular, we have not taken into ac- 
count possible correlations among the parameter uncertainties. For 
example, the uncertainties of the epoch and the longitude of the pe- 
riastron (dTo and dW) are likely correlated, and the correlation is 
probably regulated by the orbit eccentricity. 



7.1 The case of LS 5039 

The most recent est imate of the orbital parameters of LS 5039 
dAragona etaTB oQQ) are reported in the first column of Table [4] 
With these uncertainties, we find that the detection of a pulsar with 
a period faster than some seconds is, regrettably, still very unlikely. 
A search for pulsations implies sampling the space of the orbital 
parameters in order to scan the range given by the current uncer- 
tainties. The results found in this work are useful to define the max- 
imum step in the parameter sampling that guarantees a loss factor 
e 2 greater than a pre-defined value e^n- The maximum step for 
each parameter is defined as twice the maximum error calculated 
by Eqs. <52), (|55j, {64j>, (|68l>. GD, for A, W, e, T , Porb, respec- 
tively. Indeed, the maximum error dp is defined as the distance from 
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Figure 15. Power spectra of the simulated time series, demodulated varying the orbital period by dP olb = 7.6 X 10 — 4 days. The expected value of e 2 = 0.9 
is marked by the horizontal dashed line. The Coherent Power at the denominator of the y axis is the power calculated for an unperturbed demodulation 
(dP or b = 0) at the pulsar frequency u>o ■ 



Table 3. Derived formulae to calculate the maximum value of the error on the parameters in order to maintain the loss factor e 2 larger than a certain amount. 
For the maximum error de (dTo) the functions Ampl(e) and U (e) are defined in Eqs. i63\ (Eqs. (67), respectively). For the maximum error on P or b, we use 
n = T w / P or b- m the nrst column [e 2 ] 1 means the inverse function of e 2 given in the second column. For the cases we do not found an explicit formula 
for e (labeled with — ) the inverse value can be deduced by the figure listed in the third column. As explained in Section 7.1, in sampling the space of the 
parameters, the maximum step is defined as twice the maximum error given in the first column. 

Maximum error e 2 function e 2 plot 

2 (ttz) 2 ] Figure[3] 
2 (ttx) 2 ] Figure[3] 

Figure [8] 
Figure [T2l 
Figure [T4l 



dA< ^EL £ 2 = [£!^£)1 2 h. 

uioVl-e 2 cos(W) 2 L nx J L 

dw< , ±H e 2 = r £ M^)l 2 [i 

de ^[e 2 ]- 1 .{C/( e )+Ampl(e).|co B W|} 

UJqA 

ir[e 2 ]- 1 .{?7( e )- Am P' (e) ■(l+cos(2ff))} 



dT < 



^ n or b A/(i+?n~ 
10 



the true value of the parameter, where the worst case that can hap- 
pen with the sampling is that the true value of a parameter is right 
at the center of a step. Thus, the full length of the step is equal to 
2dp, and at both edges of the step, the loss factor is equal to the 
pre-defined value e m i n . In all the other cases one of the two edges 
of the step measure a loss factor e 2 > e^iir 

Table [4] shows the maximum steps that should be taken to 
search for a young pulsar with slow (300 ms), and fast (30 ms) 
period, as well as for a millisecond pulsar (with period 3 ms). The 
steps in the table are calculated such that the loss factors for the 
parameters A, e, W , and To are e p > 0.8. We will suppose to 
perform a FFT analysis with a time windows equal to the orbital 
period (T w — P or b). Therefore, the maximum steps for P or b are 
calculated setting n — 1 in Eq. {72}, and such that its loss factor is 
£p orb > 0.9. Indeed, we can be more tight with the orbital period, 
because commonly it is the parameter with the smallest uncertainty 
with respect to the others. With the steps defined in this way, using 
Eq. l(73) we see that the total loss factor is e 2 > 0.36. 

The last row of Table|4]gives the total number of trials needed 
to cover the ±1<t uncertainty ranges on the orbital parameters. For 
instance, for P psr = 300 ms the uncertainty range is covered by 
9 steps of To, 5 steps of A, 11 steps of W, and 2 steps of e. No 
trials are needed for P or b, because its ±la uncertainty ranges is 
shorter than its maximum step. The total number of trials of orbital 
parameters is then the multiplication of the former steps, N op = 



990. For all the parameters, the steps have an inverse dependence 
with the pulsar frequency. 

7.2 The case of LS I +61 303 

We calculated the maximum steps of the orbital parameters also for 
LS I +61 303, as reported in Table|5] For this system the steps of 
W, and e are much smaller than for LS 5039, because the orbit of 
LS I +61 303 is much larger. We calculated here too the possibility 
to detect a very slo w pulsar with period greater than 2 s (see, e.g., 
iTorres et al.U2012h ). 

7.3 Simulations 

To have a feeling on how the calculated steps of the orbital pa- 
rameters can allow (or not) the plausible detection of pulsations 
in gamma-ray energies, we simulated LS I +61 303 as it was ob- 
served by Perrai-LAT for 2 years, assuming a 300 ms pulsar as the 
compact object. We suppose to analyse the LAT data with an aper- 
ture photometry technique -collecting all the events w ithin an an- 
gular se paration- of 2.4 deg from the source, as done in lAbdo et al.l 
j2009bl) . This implies an event rate of 30 counts per day, and the 
background level due to the diffuse gamma-ray emission and the 
nearby sources approximately equal to 2/3 of the total counts. We 
also conservatively assume that the pulse fraction is 50%. The 
barycentered arrival time series was simulated following the same 
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Table 4. Maximum steps of the sampling of the orbital parameters calculated for LS 5039. In the last row the total number of steps needed to cover the ±l<r 
uncertainty ranges are given. 

LS 5039 



Most updated parameters Maximum steps to have the total factor e 2 > 0.36 

(Aragona eyu\^009J P pS r = 300 ms P psr = 30 ms P psr = 3 ms 



Porb (d) 


3.90608 ± 0.0001 






l.lxlO" 4 


T (MJD) 


52825.985 ± 0.053 


1.2X10- 2 


1.2xl0~ 3 


1.2xl0~ 4 


A (lt-s) 


3.33 ±0.15 


6.1xl0~ 2 


6.1xl0~ 3 


6.1xl0~ 4 


W (deg) 


236.0 ± 5.8 


1.1 


l.lxlO" 1 


l.lxlO" 2 


e 


0.337 ± 0.036 


3.5xl0~ 2 


3.5xl0~ 3 


3.5xl0~ 4 


Total number of trials 


990 


9.89xl0 6 


9.73 xlO 10 



Table 5. Maximum steps of the sampling of the orbital parameters calculated for LS I 61 +303. (The uncertainty of the orbital period is from lGregorvl fe002t) ). 
In the last row the total number of steps needed to cover the ±lcr uncertainty ranges are given. 

LS I +61 303 



Most updated parameters Maximum steps to have the total factor e 2 > 0.36 

(Aragona eyu\^009_) Ppsr > 2 s Ppsr = 300 ms P pS r = 30 ms Ppsr = 3 ms 



Porb (d) 


26.4960 ± 0.0028 






0.0013 


1.3xl0~ 4 


To (MJD) 


51057.39 ±0.23 


0.090 


1.3xl0~ 2 


1.3xl0~ 3 


1.3xl0" 4 


A (lt-s) 


20.0 ± 1.2 


0.44 


6.6xl0" 2 


6.6xl0~ 3 


6.6xl0~ 4 


W (deg) 


40.5 ± 5.7 


1.2 


1.8X10- 1 


1.8xl0~ 2 


1.8xl0" 3 


e 


0.537 ± 0.034 


0.032 


4.8xl0" 3 


4.8xi0" 4 


4.8xl0~ 5 


Total number 


of trials 


450 


1.06xl0 6 


4.41 xlO 10 


4.86xl0 15 



procedure described in Section 6.1.1, with the only difference that 
instead using Eq. d53 b for the distribution of the emitted times in 
the pulsar reference frame, we used 

P t (t) = l + »7-sm(u;o*), (74) 

where the factor 77 take into account both the background level, 
and the pulsed fraction of the source, so that in this case we have 
77 = (1 — 2/3) x 50% = 1/6. The pulsar frequency is ujo = 3.33 
Hz. The simulated arrival time series have been analysed using the 
PRESTO software. In order to demodulate the time series, and to 
fold the photons, this software need the orbital parameters as input, 
as well as an initial guess of the pulsar frequency and its derivative. 

Figure flolshows the typical plots produced by PRESTO. They 
show the most probable value of the pulsar period and its first 
derivative searched over a grid. In panel a) of Figure [T6] the de- 
modulation of the arrival time series is performed using the correct 
values of the or bital parameters, whi le in panel b) the current un- 
certainties from lAragona et alj < l2009r) are added. Whereas the first 
panel shows that the signal is clearly visible after a correct demod- 
ulation, the second regrettably shows that it is completely lost with 
the current level of uncertainty. In panel c) the demodulation is per- 
formed using the orbital parameters A, W, and e that deviate from 
the correct ones just by half of the maximum steps defined in Table 
[5] or equivalently by their maximum errors (dp max = step max /2). 
For simplicity we kept To at its true value, while the uncertainty on 
Porb is negligible in this case, as showed in Table [5] We enlarged 
the error having the maximum errors (dp max ) as units, and show 
the results we obtain in the case they deviate 2 (in panel d) and 3 
times (in panel e) the maximum errors. We can notice that the sig- 
nal is still significant (~ 6a) in panel d), even if the total factor is 
e 2 ~ 0.1. This depends on the fact that the signal to noise ratio in 
the power spectrum is proportional to the total number of counts 
observed. Specifically, P(ujo)/ P(uj 7^ wq) = Nqe 2 /2. The two 
years-long simulated observation makes possible a detection at low 



e 2 , even if the data rate and the total pulsed fraction are relatively 
small. The plots in Figure [T6] (panels c), d), and e)) show that de- 
tection is possible only if the deviations of the orbital parameters 
are within the bell shape in Figure[3]and Figure[8] depending on the 
overall pulse fraction (factor 77 in Eq.l74t. In contrast, outside that 
bell, the pulsation detection is washed out (as in panel b of Figure 

ED. 

7.4 Comparing the maximum steps with current 
uncertainties 

The comparison of the maximum steps in Tables [4] and [5] with the 
current uncertainties of the orbital parameters points out that a fine 
sampling is needed in order to detect pulsations with long observa- 
tions. In the case of a magnetar-like period, we can still hope that an 
improvement of the uncertainties with further optical observations 
can allow a detection, but for young or millisecond pulsars, sam- 
pling is unavoidable. For instance, for pulsars with a period of 300 
ms orbiting LS I +61 303, the uncertainty on the semi-major axis 
has to be smaller than 0.03 lt-s, that is ~ 40 times smaller than the 
current one. The larger the pulsar period the closer current optical 
observations are with respect to the maximum steps. The situation 
is slightly more optimistic for the consideration of a young pulsar 
in LS 5039, which has a semi-major axis of 3.33 ± 0.15 lt-s. We 
calculated that in order to have e 2 > 0.36 for a pulsar with a period 
of 300 ms the uncertainties should be improved by a factor of 10 
for To and W, a factor of 5 for A, and a factor of 2 for e. 

The study presented here is conducted assuming a pure si- 
nusoidal pulse shape for the pulsar. In contrast, the gamma-ray 
profiles show com monly two narrow peaks, or a single broad one 
(lAbdoetal.ll2010h . Therefore, the total coherent power is shared in 
several harmonics. In the periodicity search, the significance of the 
signal can be improved summing the power of different harmon- 
ics. The extension of the results in this work to higher harmonic 
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is immediate, except for the fact that in the sampling of the orbital 
parameters one should optimize the maximum steps with respect to 
the highest harmonic to sum. 

In conclusion, the analytical study presented in this paper clar- 
ifies a number of issues regarding the possible detection of pulsa- 
tions from binary systems. It is likely the basis for the development 
of a software to perform a blind search for pulsations in binary 
systems, applying an optimized sampling of the orbital parameters 
searching in different frequency ranges. 
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APPENDIX A: 

We here briefly introduce a method to evaluate the expectation 
value of the power spectrum at the signal and nearby frequencies, 

based on a statistical approach. 

Following the definition bv lScargld JT982T) . the power spec- 
trum of a sample data set {X (ti), i — 1,2, No} calculated at a 
frequency cj is given by 



P(u) 



1 

No~ 



/N \ 2 / N 



sin(u>ti 



(Al) 

Here, we consider the series of the arrival times of single events on 
a detector, i.e. the sample data set X is such that X * = 1 for each 
i. Attributing a phase value to each event (?, = ujti we have for the 
power spectrum 

2 



P ^ = No 



N 

cos6 



N 

sin#i 

i—l 



(A2) 



For a large number of events (No > 100) in Eq. l lA2t . the sums of 
the trigonometric functions of the phases, cos (Si) and sin(#i), are 
well approximated by their mean values times No. Thus, 



N 

cosOi 

(-1 



No <cos0i) = N J 



cos9- P cos (9)dcos9, (A3) 



where P cos (9) is the distribution of the values of cos# expressed as 
function of 9, and 



No (sin0i) = No / sin0-P sin (0)cfein0, (A4) 



N 

E sin (^) 

i=l 

where P s i n (8) is the distribution of the values of sin#. The expec- 
tation value of the power spectrum is obtained substituting these in 
Eq. (TAD 



P(u) = No [{cose,} 2 + (shift) 2 ] 



(A5) 



The two distributions P cos (9), and P s i n (9) can be expressed as 
function of the probability density distribution (pdf) of the phases 
attributed to each time stamp of the time series Pg{9). We use that 
the pdf of a variable z = f(x) that is function of a random vari- 
able x whose pdf P x (x) is known can be calculated as (see e.g. 
iRotondi et ai1l2004l) 



Px(x 2 j) Px(xi.i) 



f'(X2j) f'(xij) 



(A6) 
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Figure 16. Detection of the simulated pulsar in LS I +61 303. A 300 ms pulsar with 50% pulse fraction and 2 years of Femu'-LAT data have been assumed 
in the simulation. Panel a): Search for the most probable values of the pulsar period and its first derivative with PRESTO, after demodulating the arrival time 
series with the correct orbital parameters. Panel b): The same but after demo dulating the arrival time series using parameters that deviate from the correct 
ones by an amount equal to the current uncertainties in|Araeona etal.1 j2009h . Panel c): The same but with the demodulation performed using the orbital 
parameters A, e, and W that deviate from the correct ones by half of the maximum steps defined in table \5\ (step max ), or equivalently by their maximum 
eiTors dpmax = step max /2. Panel d): The same but the orbital parameters A, e, and W deviate by 2 times the maximum errors dp max . Panel d): and 3 times 
the maximum errors c£p ma x- 



Here, the prime (as in /') represents a derivative with respect to x, 
and the intervals [xij , X2.i] are those for which for a given z — zq, 
f(x) < zq. In our case, P z in Eq. iA6\ corresponds to P cos or 
P s i n , while P x is the phase distribution Pg. A detailed calculation 
of Pcos and P s i n leads to the result that both P cos and P s i n can be 
expressed as a linear function of the sums 



APPENDIX B: 



In Section 3, we approximated the eccentricity with E ~ r2 or b.i e 
within the perturbation functions in Eq. d 1 3 fc . In this Appendix we 
shall evaluate the level of accuracy of this approximation. 



N-l 

J2 Pe(2n(i + k)±e), (A7) 

1 = 



where k — 0, 1/2, 1, and N is the number of rotations made by 
the plausible neutron star during the whole observation T b s . Since 
in Eq. ( IA7t . ^ 9 < 2n, it corresponds to the distribution of the 
phases folded by 2ir. Substituting the solutions of P cos and P sln 
in Eqs. dA3t and JA4b respectively, and these in Eq. dA5b . we ob- 
tain the direct dependence of the power spectrum by the pdf of the 
phases Pg(8). 



The worst possible case for the approximation is when the or- 
bit has an eccentricity close to 1 . Figure IB II shows how different 
are the correct values of E(t e ) computed numerically inverting the 
Eq. d 1 3 b . with respect to its approximation. In the same figure we 
show the differences between the sine of E(t e ) and the sine of its 
approximation (hereafter E). The linear and the sine dependencies 
on E are the two instances concerned by the perturbation function. 
The former appears only in the case the perturbed function is due 
to the uncertainty on the orbital period. In the case of the linear 
dependence, the larger is Q OT ht e , the smaller is the relative dif- 
ference between E(t e ) and E. So the approximation improves for 
long observations, and the difference becomes negligible by setting 
the zero of the time stamps t e far enough from the observation time 
In the case of the sine dependence, we estimated the mean squared 
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Figure Bl. Top: Numerically computed eccentric anomaly E/2tt (in blue), 
compared with its approximation E/2ir = Q OI ^,te/2ir (in red). Bottom: 
Numerically computed sine of the eccentric anomaly (in blue), compared 
with the sine of the approximation adopted (sin E, in blue). Magnitudes in 
both panels are shown as a function of f2 or ] 3 t e /27r. 



difference using the following formula 



X = 



\ 



sini? — sini5 



(Bl) 



using a sampling N s = 1000 of f2 or bi e in the range 0, 2n. In the 
case showed in figure lBl"l (e = 1) x = 0.346. Figure lB2l plots \ 
as calculated for several values of the eccentricity. All in all, these 
approximations do not seem, a priori, to introduce a significant de- 
viation. This is later confirmed by simulations. 



APPENDIX C: 

In this Appendix we find the conditions under which the func- 
tion tew = t e — f(t e ) increases monotonically. We are going to 
demonstrate that these conditions are satisfied by all the binary sys- 
tems that are already observed. To increase monotonically, the first 
derivative of t ew with respect to t e should satisfy 



f-< L 

at e 

In Eq. (O the perturbation function is defined as 

f(E(t e )) = g p dp, 



(CI) 



(C2) 



Figure B2. Mean squared difference x between the sine of the eccentric 
anomaly and the sine of its approximation (Eq. lBU . calculated for several 
values of the eccentricity. 



where p is the orbital parameter to which the perturbation function 
is referred to. We stress that in Eq. dC2t both / and Ar are func- 
tions of the eccentric anomaly E, which in turn is function of t e . 
Taking this into account, the first derivative of / can be calculated 

as 



df(E(t e )) = d 2 A R dp = d_ 
dt e dt e dp dp 



dAi 

dt e 



dp, 



(C3) 



where the exchange of the order of the derivatives in the last pas- 
sage is allowed because t e and the orbital parameters p are inde- 
pendent variables. Since A« is only indirectly a function of t e , its 
derivative with respect to it is 



8Ar _ dA R dE 

dt e ~ dE dte 

From Eq. Q, the derivative of A_r respect to E is equal to 

dA R 



dE 



= M(A, e, W)cos{E + <j>(e, W)), 



(C4) 



(C5) 



where here we explicitly write the dependence of M and cj> (de- 
fined in Eq.[8]l on the canonical orbital parameters. To evaluate the 
derivative of the eccentric anomaly E with respect to t e , we can 
differentiate Eq. d 1 3 b 



dE 

dt e 

from which we obtain 



e cosE 



dE_ 

dt e 



dE 

dt e 



1 — e cosE 

Substituting Eqs. tC4llC31IC7t in Eq. O, we get 

df(E(t e )) _ 



(C6) 



(C7) 



dt e 



d_ 

dp 



e cosE 



M(A, e, W)cos{E + (f>(e, W)) 



(C8) 



We will find useful to indicate the expression within the square 
brackets in Eq. |C8| as the function = g{p, E). To calculate the 
derivative in Eq. dC8t , and to evaluate the inequality df/dt e < 1, 
we need to analyze case by case each orbital parameter p. 
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CI [p = A] 

For the projection of the semi-major axis (A), Eq. iCSl becomes 

df{E{t e )) _ 



JL 

dA 



dt e 
^orb 

1 — e cosP 

^orb 



1 — e cosP 



M(A,e,W)cos(E + <j>(e,W)) dA = 

cos(E + <f>)dA < 1, (C9) 



yl — e 2 cos 2 H / 



where the last inequality is what we need to prove. We eliminate 
the dependence of the uncertainty dA by the eccentric anomaly E, 
because the former is the uncertainty on a parameter that describe 
the system, and is fixed in time. Since cos(P + (j>) has as maximum 
and minimum values 1 and —1, respectively, and the term in square 
bracket is positive, dA is still well constrained by 

fiorb 



1 — e cosP 



\/l — e 2 cos 2 W / 



\dA\ < 1. 



(CIO) 



In order to be as conservative as possible with respect the un- 
certainty dA, we can maximize the term in square brackets in 
Eq. (fCTOj. With this purpose the denominator (1 — e cosP) is min- 
imized setting cosE = 1, while the root square at the nominator is 
maximized setting cosVF = 0. We so obtain 

1 - e 



\dA\ < 



We can multiply and divide by (1 + e), obtaining 

1-e 2 



\dA\ < 



(Cll) 



(C12) 



n orb (l + e) ' 

We now minimize the right hand side, so that the constraints on the 
uncertainty are the most conservative possible, getting, 

Porb(l-e 2 ) 



\dA\ < 



4tt 



(C13) 



where we considered that Sl or b = 2n/ P or b ■ 

It may seems at this point that the inequality would not hold, 
since for the two body problem with a gravitational potential the ec- 
centricity could be as close as possible to 1. But tidal forces would 
come to help. 

Tidal forces tend to circularize the orbit of a binary system, 
affecting the period and the eccentricity at once. They strongly de- 
pend by the separation (r) of the objects in the system, varying 
as 1/r 6 dLecar. Wheeler. & McKedl 19761) . Thus, in systems with 
short periods, tidal forces are strong. On the other hand, in systems 
with long periods, they are still not negligible when high eccentrici- 
ties make the two objects very close at periastron. Systems affected 
by tidal forces modify in short time their orbit until they reach a bal- 
anced status in which either the se forces ar e no lo nger important, 
or the orbit become circular. In lHalbwachsl d2008l) . the maximum 
eccentricity that a system can have without being significantly af- 
fected by tidal forces is expressed as function of its period by 



f'Ma 



1 - (fboff/Po, 



,2/3. 



(C14) 



where Pcoff is the circularization limit, meaning that systems with 
periods that are shorter than the latter inevitably fall in a circu- 
lar orbit. Using th e sample of the observed spectroscopic binaries, 
lHalbwachsl ( f2008l) (and references therein) estimate Pcoff = 5-j-10 
days. These values are compatible w ith the sample of the HMXBs 
analyzed by iTownsend et al" (compare their figure 6 with 

Eq. lC14b . To be conservative, in this work we will consider Pcoff = 
5 days. 



The right term in Eq. dC13t is further minimized if the maxi- 
mum eccentricity expressed by Eq. dC14t is taken into account. So 
that we have 



\dA\ 
A 

\dA\ 
A 



< 



P C o« 



2/3 



Porb 

4^4' 



for Porb > -Pcoff ; 



(C15) 



< 



for P orb < Pcoff, 



(C16) 



Porb 
-Porb 

4^4' 

where we have also divided both the left and right side by A, to take 
into account the relative uncertainty |<L4|/A In order to express 
the right hand side of the inequalities as function only of the orbital 
period, we can invoke the third Kepler's law 

1/3 



A 



G 



Att 2 ^ {Ms + M ps 



(C17) 



where G is the gravitational constant, c the speed of light, M psr the 
mass of the compact object, M B the mass of the companion star, 
and sin i is the inclination of the orbit plane respect to the observer. 
When the companion star is very massive, the ratio M^/{M S + 
M psr ) 2 ~ M s , while when the masses are similar Mg /(M s + 
M psr ) 2 ~ M s /4. Since A is at the denominator of the terms in the 
right hand side of Eqs. JC 151 IC 1 6b . the most stringent constraint on 
dA can be found by setting sin i = 1, and M^ / (M s + M psr ) 2 — 
M s in all the cases. Finally we have 



M<9.54xl0 2 
A 



for P orb > Pcofr; 

M <9 . 54 xl0 2 
A 

for Porb < Pcoff- 



' Ms ' 
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" Pcoff " 
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2M Q _ 
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' Ms ' 


-1/3 


Porb 


1/3 


_2M _ 




5days 





(C18) 



(C19) 



From Eqs. JC 1 81 |C19l l, we can clearly see that for low mass sys- 
tems with orbital periods of few days, the uncertainty dA is to- 
tally unconstrained, since it is only required that the relative un- 
certainty |dA|/j4 < 10 , which is obviously true for uncertainties 
that make sense. In general, we can say that a constrain on dA is 
effectively worrisome when |dj4|/A < 1. This happens only for 
systems with very large masses (M s > 50Mq), and unrealistic 
periods as long as 10 6 years (from Eq. IC18t . or as short as few 
milliseconds (from Eq. |C19| l. In other words the monotonically- 
increasing-function condition for the function t ew — t e — f(t e ) is 
always satisfied for all the observed binary systems when the per- 
turbation function is due to the uncertainty on the projection of the 
semi-major axis. 

C2 [p = W] 

For the longitude of the periastron (W), Eq. jC81 > (with the imposi- 
tion of the searched inequality) reads 

df(E(U)) _ 



dt e 



(C20) 



d 
dW 



Q 



orb 



1 — e COSp 



M(A, e, W)cos{E + cj>(e, W)) 



dW 



0. o 



e cosP 



dM 

^L. COS (E + ^e,W))- 



M(A,e,W)^ S in(E + c/>(e,W)) 



dW < 1. 



Within the square brackets in the last equation, we can apply the 
trigonometric identity described in Eq. ® to cos(P + <f>) and 
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sin(P + 



, obtaining 



1 - ecosE \j \dW ) \dW ) 

sm(E + <f> + /3)dW < 1, 



(C21) 



where /3 is a function of A, e, and W. Since the factor multiplying 
sin(_E + (f> + /?) is positive, and the sine oscillates between 1 and 
— 1, the constraints on dW can be written as 



\dW\- 



fiort 



1 — e cos_E 

Maximizing the factors that multiply \dW\ we obtain 



< 1. 



(C22) 



n 



1 — e cosS 



< 



1 + e 2fi orb 



l 



l 



M(A, e, W) = A\/l~ e 2 cos 2 W < A, 



(C23) 



dM 
dW 



Ae 2 sinW cosW 



e 2 cos 2 W 
l + tan(W) 2 ) 



2y/T^l? 
1 



dW 1 - e 2 +tan(W) 2 

Eq. dC22t then becomes 

Porb (1 



\dW\ < 



e 2\(3/2) 



(C24) 



2tt^/5 A 

where we have substituted Q ib = 2n/ P OT b- Substituting the ec- 
centricity e by 6Max as defined in Eq. l IC 14b in Eq. lC24l we get 

Pcoff 1 



\dW\ 
\dW\ 



2-kJZA' 

Porb 1 



for P or b > Pcofi; 
for P orb < Pcoff- 



(C25) 
(C26) 



271-^5^4' 

Finally, introducing the third Kepler's law (Eq, IC17l > with sin i = 1, 



and Mi /{M a + M ps 
\dW\ < 8.54 x 10 2 



M s , we get 



M s 



2M P 



-1/3 r 



PCoff 



5days 



-i -2/3 



5days 



for P orb > Pcoff; (C27) 



|dW| < 8.54 x 10 2 



" M s ' 


-1/3 


Porb 


1/3 


_2M _ 




5days 





for P orb < Pcoff. (C28) 



Since W is an angle, the constraints on its uncertainty dW are ef- 
fective if dW < 2ir rad is required. A greater uncertainty means 
that we do not have any knowledge of W whatsoever. From the 
equations above, we can calculate that dW is effectively con- 
strained only for binary systems with periods longer than 7 years 
and masses M, > 20Mq. A more constraining but still very 
large uncertainty dW < 1 rad satisfy the monotonic condition of 
the function t ew for all the binary systems with M s < 50Mq, 
and period shorter than 70 years. Thus, it is safe to consider that 
the monotonically-increasing-function condition for the function 
tew = t e — f(t e ) is always satisfied for the systems of interest. 



C3 



e] 



Since the eccentric anomaly depends by the eccentricity (e), for this 
parameter Eq. l |C8l > (with the imposition of the searched inequality) 
reads 

d 



df(E(t e )) 

dte 



()(: 



b (e,P)He=g.l + |§f d e<l (C29) 



where we defined with g(e, E) the expression in the square brack- 
ets of Eq. JC8t . Term by term we have 



dg_ 

de 



or b 



e cosi? 



AfcosP dM 

^ ~de~ 



1 — e cosE 



cos(E - 



M-^-sm{E + 4>) \ (C30) 



Ho 



BE 



dE_ 

de 



1 — e cosE 



e M sinP 



e cosP 



cos(P + 4>)- 
Msia(E + </>)} 



e cosP 



(C31) 
(C32) 



Substituting Eqs. JC30l dC3ll and dC32l in Eq. tC29\ , and applying 
the trigonometric identity described in Eq. {6]l to cos(i5 + (j>) and 
sin(P + <jj), we obtain 



n 



orb 



1 — e cosiS 



McosE + dM eMsin 2 E 
1 — ecosP de (1 — ecosE) 2 



M d(j> MsinE 
de 1 — e cosP 



sin(P + <j> + p)de < 1 (C33) 



where j3 is a function of A, e, and W. The terms in the brace brack- 
ets can be maximised using the following inequalities, in addition 
to those in Eq. dC23t 



( dM\ 2 _ 
[de) 



Ae cos W 



VI -e 2 cos 2 W 
etan(W) 



Ae 



VT 

< 



(f-e 2 +tan(P7) 2 )VT^ 2 2(1 - e 2 ) 



The constraint on de is then reduced to 

Porb (1 - 



\de\ < 



27rV61 



A 



(C34) 



(C35) 



(C36) 



where we have already substituted f2 rb = 27r/P or b. Taking into 
account Eq. JC14b . and the third Kepler's law (Eq. lC17t . we finally 
get 
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for P orb > Pcoft; 



\de\ < 2.44 x 10 2 
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(C37) 



(C38) 



Since the eccentricity is in the range ^ e < 1, its uncertainty 
de is really constrained if de < 1. Taking this into account, we 
can calculate from the equations above that also in this case the 
monotonically-increasing-function condition for t ew = t e — f(t s ) 
is satisfied for all realistic binary systems. 



C4 [p = T ] 

Since the eccentric anomaly depends by the epoch of the perias- 
tron (To), for this parameter Eq. dC8t (with the imposition of the 
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searched inequality) reads 



(C39) 



where we defined with g (To, E) the expression in the square brack- 
ets of Eq. JC8t . Term by term we have 



dg 

on 



= 



(C40) 



"ort 



(C41) 



dE _ 

dfo ~ 1 - e cosE 

while is given in Eq. JC3U . Substituting in Eq. dC39t . and ap- 
plying the trigonometric identity described in Eq. ® to cos(E-\-<f>) 
and sin(£' + </>), we obtain 



(1 - ecosP) 2 



e M sinP 



1 — e cos i5 

sin(P + 4> + P)dT < 1 



(C42) 



where /? is a function of A, e, and W. Maximising the terms in the 
brace brackets using the inequalities in Eq. dC23t the constraint on 
de is then reduced to 

1-e 2 ) 3 



I dr 1 

-Porb 



< 



Porh 

16tt 2 VE A 



(C43) 



where we have already substituted Q ib = 2n/P OI b- Taking into 
account Eq. JC 14b . and the third Kepler's law (Eq, IC17l >. we finally 
get 
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A constrain on dTo is effective when \dTo\/ P OY b < 1. This hap- 
pens only for systems with very large masses, and unrealistic very 
long periods. In conclusion, the monotonically-increasing-f unction 
condition t ew = t e — f(t e ) is satisfied. 

C5 [p = P orb ] 

Since the eccentric anomaly depends by the orbital period (P or b), 
for this parameter Eq. JC81 > (with the imposition of the searched 
inequality) reads 



df(E(t e )) 



d 



dt e 
dg 



dD. OI 



1 + 



dg dE 



[g(P OTh ,E)] dP orh 



dE dn o 



dp m 



-dPorb < 1 



(C46) 



where we defined with g(Q. OY b,E) the expression in the square 
brackets of Eq. dC8t , and we take into account that Q. ora = 
2n I P or b Term by term we have 

d^lorb 



2-k 

p\ 

orb 



dPorb 

dg _ Mcos(E + <f> 
dflorb 1 — e cosP 



(C47) 



(C48) 



dE 

dfl OT y_ 



To 



e cosP 



(C49) 



while ^ is given in Eq. jC3U . Substituting in Eq. dC46fc . and ap- 
plying the trigonometric identity described in Eq. ((6) to cos(P + </>) 
and sin(P + </>), we obtain 



2nM 



(1 - ecosP)P Q] 



1 - 



f2 or b esinP (t e — T ) 



+ 



fiorb (te — To) 



(1 - ecosP) 2 

sin(P + 4> + P) < 1 (C50) 



e cosP 



where /3 is a function of A, e, and W. Maximising the terms in the 
brace brackets using the inequalities in Eq. dC23b the constraint on 
dPorb is then reduced to 



IdP,, 



(1-e 2 



D orb 167T\/l + 47r 2 (t e -T ) 2 /P c 



orb 



A 



(C51) 



where we have already substituted f2 rb = 2n/ P or b . If To is set 
in order to be as close as possible to the beginning of the observa- 
tion or better within it, then (t e — To) ^ T b s . The inequality in 
Eq. dC51l > is still valid if we substitute (t e — To) with T b s = raP or b- 
Taking into account Eq. JC 14b . and the third Kepler's law (Eq. 
|C17t , we finally get 



|dPorb| 2.4 x 10 2 



Porb VI + 47T 2 n 2 



Ms 



2Mr, 



-1/3 



PCoff 



5days 
for P orb > Pcoff 



5days 



-5/3 



(C52) 
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(C53) 



Commonly, the uncertainty on P or b is such that |dP rb|/Porb < 
(C44) 10 -3 . This level of constrain is reached in Eq. dC52| l and l |C53t only 
for n > 10 4 . On the other hand, in the analysis proposed use time 
windows T w — P or b for long period systems (n — 1), or for orbital 
periods of few hours a time window of a week implies n to be of 
(C45) the order of 100. A part for n, is not satisfied only for systems with 
very large masses, and unrealistic very long periods. In conclusion, 
the monotonically-increasing-function condition t ew = t e — f(t s ) 
is satisfied. 



C6 The most general case 

So far we have considered the case in which only a single param- 
eter is estimated with a given uncertainty, while all the others are 
known, so that their contribution is null to the non-factorised per- 
turbation function, as defined by Eq. © 



/ 



V- dA Rj 



(C54) 



In practice, however, all orbital parameters have an uncertainty, and 
the condition over the function t W e must be satisfied when the full 
derivative of Eq. \C54l with respect to i e (that is equal to the sum 
of the terms in Eq. IC3I relative to each one of the parameters) is 
lower than 1, i.e., 



df(E(t e )) 

dt e 



d fdA R 

^ d P \ dte 



dp < 1. 



(C55) 



This inequality is satisfied if we divide by 5 each of the formerly 
derived constraints, i.e., Eqs. ( tCl^[Cl^[C27llC28IIC37|[C3^1IC44l 
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IC45IIC52llC53l >. This introduce no significant worry for any of the 
parameters, since the conditions were easily satisfied in all cases of 
realistic binary systems. 
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